cc Contains the subroutines of the piecewise perturbation methods
cc to solve the radial Schroedinger eq. for complex wavenumber.
cc Written by
cc                   L. Gr. Ixaru, M. Rizea,
cc  Institute of Physics and Nuclear Engineering "Horia Hulubei"
cc                            Bucharest
cc                               and 
cc                            T. Vertse
cc Institute of Nuclear Research of the Hungarian Academy of Sciences
cc                             Debrecen
	subroutine swf(rrxmax,rixass,rrstep,
     1                 rmatch,mnbo,icore,
     2                 npar,at,zt,v_0,r_0,a,v_so,
     3                 rolrtol,
     4                 igamm,
     5                 lbas,jbas,esw,
     6                 scawf,
     7	               ck,rgaussx1,u,up,ngauss,v_ws)
	parameter(nmax=2000,maxs=250,neqd=nmax,mni=3000,mnia=300)
	parameter(mxrd=1000,maxd=500,nkmax=30,nbase=1000,maxc=nbase)
	implicit double complex (a-h,o-q,s-z)
	implicit double precision (r)
cc040415
	real*8 at,zt,v_0,r_0,a,v_so
	dimension v_0(2),v_so(2)
cc040415
	dimension nzp(2)
	common/sabsa/rtol,rmax,rstep,nstep,ll,j2,nzp,nkval
c	common/rolrtol/ rolrtol
c	common/acsw/scawf(mxrd)
cc040415
	dimension scawf(mxrd)
cc for k^2 cal. cubic sample line fit coeff. 	
	dimension b(neqd),c(neqd),d(neqd)
	dimension up(1:*),u(1:*),rgaussx1(1:*),v_ws(1:*)
cc	
	common /spcalc/ ispn,isps
	common/poti1/v1s,v10,v11,v12
	common/poti2/pas(mni),v(mni),v1(mni),v2(mni),intf
	common/poti3/apas(mnia),av(mnia),av1(mnia),av2(mnia),intaf
	common/genwork/par(3),rl
	common/radial/ renfac,retak,lorb,jorb
	common /abscis/ rxeq(0:nmax),rx(0:neqd),nbo,numeq,numneq
	common/xxxx/rxx,rpsimin,rffi,rbestifi
	common/intdia/NUMB,MM
	common/radstep/rhs
cc040415
c	common/c1pot/rrxmax,rixass,rrstep
c	common/c2pot/rmatch,mnbo,icore
c	common/c3pot/npar(2),at,zt,v_0(2),r_0,a,v_so(2)
	common/exc3pot/attt,ztt,v_00(2),r_00,aa,v_soo(2)
	real*8 attt,ztt,v_00,r_00,aa,v_soo
cc040503
	common/rolo/igam
cc040503
cc040415
	dimension y(0:nmax),yp(0:nmax)
	dimension sy(0:nmax),syp(0:nmax),dy1(0:nmax),dy2(0:nmax),
     & dyp2(0:nmax),dyp1(0:nmax)
	external spot,pot,apot,aspot
	data czero/(0.d0,0.d0)/,cunit/(1.d0,0.d0)/
c	write(*,*)' l,j,e=',lbas,jbas,esw
cc040503
	igam=igamm
cc040503
cc040415
	attt=at
	ztt=zt
        v_00(1)=v_0(1)
        v_00(2)=v_0(2)
        r_00=r_0
        aa=a
        v_soo(1)=v_so(1)
        v_soo(2)=v_so(2)
cc040415
	rixmax=0.d0
	rrxass=rrxmax
	rstep=rrstep
	nbo=mnbo
	rxx=0.1
	par(2)=dcmplx(rrxmax,rixmax)
	par(3)=dcmplx(rrxass,rixass)
	rhs=rstep
c       write(*,*)' par=',par(2),par(3)
cc numeq = number of equidistant mesh points
cc rxeq(i)  = equidistant mesh points
    	numeq=rrxmax/rstep
	if(numeq.gt.mxrd)then
	 write(2,*)' numeq.gt.mxrd, increase rstep'
	 stop
	endif
	rstep=rrxmax/numeq
	nstep=numeq
	rmax=rrxmax
	rrxmin=nbo*rstep
	rixmin=0.d0
	par(1)=dcmplx(rrxmin,rixmin)
	rxeq(0)=0.d0
	do int=1,numeq
	 rxeq(int)=int*rstep
c	write(*,*)' rxeq=',rxeq(int)
	enddo
crolo	rtol=1.d-10
	rtol=rolrtol
	reps=0.001d0
	if(abs(reps).lt.1.e-3)reps=0.1d0        
	lsel=lbas
	jsel=jbas
	isps=0
	ispn=0
cc the information of whether the potentials are furnished as functions
cc or in tabulated form is settled        
	lorb=lsel
	jorb=jsel
	rl=lorb
c	rxx=rstep
c       rxx=0.1
	nzp(icore)=icore
cc generating potential terms
	itz=icore
	call potgen(lbas,itz)
c       at this moment all information on the potentials is available 
c       and the numerical solution can start.
c       first,the ends of the intervals i1,i2 and the value of the assymptotic
c       point must be settled. they are stored in par(i),i=1,2,3,resp.         
c        
c       the program now generates the the weighted potential values v1s,v10,
c       v11 and v12 which are further necessary for the integration on i1.
c  	write(*,*)' before spi1'
	call spi1(par(1),v1s,v10,v11,v12)
c   	write(*,*)' after spi1'       
c       data which are further necessary for the integration on i2,viz. 
c       partition of i2 (vector pas), the three weighted potential values
c       (vectors v,v1 and v2) and the total number of the intervals intf 
c       are now generated. the user should provide the maximal value allowed
c       for intf. this governs the dimension of the vectors pas, v, v1 and v2.
c       in output, intf contains the number of intervals which are actually
c       used. 
        intf=mni
c	  	write(*,*)' before sli'
        call sli(rtol,pas,v,v1,v2,par,rl,spot,pot,intf)
c	  	write(*,*)' before spi1'
c       
c       the same data are now constructed to cover the ray between par(2)
c       and par(3), i.e. in the assymptotic region. they are stored in vectors
c       apas,av,av1 and av2 and in integer intaf.
c        
        intaf=mnia
        call sli(rtol,apas,av,av1,av2,par(2),rl,aspot,apot,intaf)
c       just an illustration of how many intervals are used on i2 and on i3
c       follows        

c        
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c  numneq = number of non-equidistant mesh-points
c  rx(i),i=0,1,..numneq  = non-equidistant mesh-points
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	numneq=intf+nbo
	rx(0)=0.d0
	do int=1,nbo
	rx(int)=int*rstep
	enddo
        x=par(1)
        do 16 int=1,intf
        x=x+pas(int)
	rx(int+nbo)=x
16      continue
50      format(1x,'int= ',i5,1x,' x= ',2f8.5,' pas= ',2f8.5,2x,2d10.3)

c       the value of the matching interval(i.e. the matching point will be
c       located at  the r.h.end of the interval labelled by match in i2 for
c       all calculations which follow) is settled; in this illustrative
c       example it is settled as the interval which contains the point x=7.5.
c
        x=par(1)
        do 8 int=1,intf
        x=x+pas(int)
        if(abs(x).gt.rmatch) go to 9
8       continue
9       match=int
        newit=10

	do i=1,ngauss
	   rxtemp=rgaussx1(i)
	   v_ws(i)=pot(rxtemp)+spot(rxtemp)/rxtemp
	end do
cc040503
c	if(igam.eq.0)goto 3923
	if(igam.eq.1)goto 3923
cc040503
c
c       the procedure to locate the eigenvalues follows. the user should first
c       give a first guess for the complex k  or the energy esp
       
c  loop for the starting values for the discrete basis states
	energy=esw
	ck2=renfac*energy
cc040503
c        ck=cwn(ck2)  
        ck=cwn(igam,ck2)  
cc040503
c        write(*,*)' e before',energy
	call sfindev(ck,match,swght,gap,dwght,rtol,ierr)
	
	if(ierr.ne.0) then
	print *,ierr
	print *,' ierr.ne.0, I stop'
	stop
	endif
c        k eigenvalue is found and stored in ccwk
	energy=ck*ck/renfac
c        write(*,*)' e after',energy
c       we now have the eigenvalue in the output ck, and also the values of
c       three parameters (swght, gap, gwght) which are further used in the 
c       computation of the normalized eigenfunction corresponding to this 
c       eigenvalue.
c        
c
c       The calculation of the normalized eigenfunction follows.
c       First the eigenfunction is calculated at the original, nonequidistant
c       mesh. Later it is interpolated and written out at equidistant
c       mesh points.
c
c       nbo is the number of equidistant mesh points in the interval I1.
c       The values of the eigenfunctions and its derivative wrt. r
c       at points rx(i), i=1,2,..,nbo,...,numneq  are stored
c       in vectors y  and yp respectively.
c
        call snorm(nbo,ck,gap,match,swght,dwght,y,yp)
c       write out at the original mesh	
c       interpolate and write out at equidistant mesh
        call eqmesh(y(nbo),dy1(1),yp(nbo),yp(numneq),dy2)
c        write(*,*)' normalised solution follows'
c rolorolo
c        write(12,871)energy,ck
c871     format(2f10.3,2f20.13)
        do ir=1,numeq
	if(ir.le.nbo) then
c	ubase(ir,ibasis)=y(ir)
	scawf(ir)=y(ir)
	else
c	ubase(ir,ibasis)=dy2(ir)
	scawf(ir)=dy2(ir)
	endif
	enddo
c        do ir=2,numeq,2
c rolorolo
c        write(12,834)rstep*ir,ubase(ir,ibasis)
c834     format(f10.3,2d22.14)
c        enddo
	esw=energy
cc store yp in up for k^2 cal.
	coul=2.d0*retak
	call spline(rx(1:numneq),yp(1:numneq),b,c,d,numneq)
	do i=1,ngauss
	   up(i)=cispline(rgaussx1(i),rx(1:numneq),yp(1:numneq)
     &	   ,b,c,d,numneq)
	end do
	
	call spline(rx(1:numneq),y(1:numneq),b,c,d,numneq)
	do i=1,ngauss
	   u(i)=cispline(rgaussx1(i),rx(1:numneq),y(1:numneq),
     &   b,c,d,numneq)
	end do
	
	return
3923	continue
c
c       calculation of the scattering states 
c 
c  ffi is the rotation angle of r for bound and resonant states
c  and also for the outgoing components of the scattering states
        rffi=datan(1.d0)-0.5d0*rpsimin
c  bestifi is the rotation angle for  the incoming components of 
c  the scattering states
        rbestifi=-2.d0*datan(1.d0)+rpsimin
c calculates scattering solutions along the contour
c	call scatt(match,lsel,jsel,nzp,rstep,rmax,esw)
	call scatt(match,lsel,jsel,nzp,rstep,rmax,esw,scawf,
     &	          ck,rgaussx1,u,up,ngauss)
	itz=icore
        nzp(icore)=itz
	return
      end
c===============SUBROUTINE SCATT======================
c	subroutine scatt(match,lsel,jsel,nzp,rstep,rmax,esw)
	subroutine scatt(match,lsel,jsel,nzp,rstep,rmax,esw,scawf,
     &	ck,rgaussx1,u,up,ngauss)
	implicit double precision(r)
	implicit double complex (a-h,o-q,s-z)
	parameter(maxs=250,nseg=20,nmax=2000,neqd=nmax,nbase=1000)
c        parameter(mxrd=150,mni=1500,mnia=300)
        parameter(mxrd=1000,mni=3000,mnia=300)
	double precision alfa,aaaa,bbbb,cccc,dddd,absc(100),
     & wei(100)
        dimension nzp(2)
	dimension sy(0:nmax),syp(0:nmax),dy1(0:nmax),dy2(0:nmax),
     & dyp2(0:nmax),dyp1(0:nmax)
	double complex delt,ddelt,delta,delta0,deltaf,test(0:20,2)
     &             ,enec,emin,emax,zi,zf
cc
c        common/acsw/scawf(mxrd)
        dimension scawf(mxrd)
cc

cc for k^2 cal. cubic sample line fit coeff. 	
	dimension b(neqd),c(neqd),d(neqd)
	dimension up(1:*),u(1:*),rgaussx1(1:*)
	
        common/genwork/par(3),rl
        common/tempwf/sy,syp,dy1,dy2,dyp1,dyp2
	common/radial/ renfac,retak,lorb,jorb
        common/poti2/pas(mni),v(mni),v1(mni),v2(mni),intf
cc040503
	common/rolo/igam
cc040503        
	common/zcont/zi(nseg,0:10,2),deltae(nseg,0:10,2),
	1 alfa(nseg,0:10,2),stepc(nseg,0:10,2),
	2 ndivc(nseg,0:10,2),mint(0:10,2)
        common/scatter/cck(maxs),us(nmax,maxs)
	common /abscis/ rxeq(0:nmax),rx(0:neqd),nbo,numeq,numneq
	common/radstep/rhs

        l=lsel
        jj=jsel
	indj=1
	if(jj.lt.2*l)indj=2
c
        pi=4.d0*atan(1.d0)
	enec=esw
        ck2=renfac*enec
cc040503        
c        ck=cwn(ck2)
        ck=cwn(igam,ck2)
cc040503        
c        ck=enec
c        enec=ck**2/renfac
c
c       complex k is in ck  and calculation of the normalized scattering
c       solution by calling subroutine scatwf as follows. 
c       you will have the wf and its
c       derivative in vectors sy and syp. 
c       The vectors dy1, dyp1, dy2 and dyp2
c       play here the role of working vectors.
c
        call scatwf(nbo,ck,match,sy,syp,dy1,dyp1,dy2,dyp2)
	rx(0)=0.d0
	do int=1,nbo
	rx(int)=int*rstep
	enddo
        x=par(1)
        do 16 int=1,intf
        x=x+pas(int)
	rx(int+nbo)=x
16      continue
c
c
c       now the wavefunction at equidistant partition is generated.
c
        call eqmesh(sy(nbo),dy1(1),syp(nbo),syp(numneq),dy2)
        do ir=1,numeq
	if(ir.le.nbo) then
	scawf(ir)=sy(ir)
	else
	scawf(ir)=dy2(ir)
	endif
	enddo
c
c  set correct phase of the radial wfn
c
ccc	print *,' ck',ck
        cphase=(1.d0,0.d0)
        if(dreal(scawf(1)).lt.0.d0)cphase=(-1.d0,0.d0)
        do k=1,numeq
        scawf(k)=scawf(k)*cphase
        enddo
51	format(2x,f7.4,2x,d20.10,2x,d20.10)
91	format(2x,f7.4,2x,d20.10)

cc for k^2 cal.
        call spline(rx(1:numneq),syp(1:numneq),b,c,d,numneq)
	do i=1,ngauss
	   up(i)=cispline(rgaussx1(i),rx(1:numneq),syp(1:numneq),
     &   b,c,d,numneq)
            up(i)=up(i) *cphase
	end do
	
	call spline(rx(1:numneq),sy(1:numneq),b,c,d,numneq)
	do i=1,ngauss
	   u(i)=cispline(rgaussx1(i),rx(1:numneq),sy(1:numneq)
     &	   ,b,c,d,numneq)
            u(i)=u(i)*cphase
	end do

	return
	end
cc040503        
c============FUNCTION CWN====================
	complex*16 function cwn(i,e)
	implicit real*8 (a-h,o-z)
c
c	this function calculates the complex wavenumber cwn from the
c	square of the complex energy  times 2m/h**2 (k**2=e*2m/h**2)
	complex*16 e
c	print*,' igam in cwn',i
c	print*,' energy in cwn',e
	cwn=sqrt(e)
c	print*,' k in cwn before if',cwn
	if(dimag(cwn).eq.0)return
	if(dimag(cwn).gt.0)then
	if(i.eq.-1)return
	cwn=-cwn
	endif
	if(i.eq.-1)cwn=-cwn
c	print*,' k in cwn after if',cwn
	return
	end	
cc040503        
cc040503        
c============ old FUNCTION CWN====================
c	complex*16 function cwn(e)
c	implicit real*8 (a-h,o-z)
cc
cc	this function calculates the complex wavenumber cwn from the
cc	square of the complex wavenumber e being proportional to
cc	the complex energy.
cc	the complex e is taken from a combination of the physical and
cc	unphysical riemann sheets:
cc	from the physical sheet the upper half including the negative
cc	energy axis is taken (so the bound states are included)
cc	from the unphysical sheet the half which was not considered on
cc	the physical sheet. (so decaying resonances are included but
cc	capturing ones and antibound states are not included)
cc
cc	note: on vax csqrt(e) would give the same value except in the case of
cc	bound states where it might give the root corresponding to the anti-
cc	bound state if the e has a small negative imaginary part.
cc
c	complex*16 e
c	common/xxxx/xx,psimin,ffi,bestifi
c	data twopi/6.283185307179586d0/ pi/3.14159265358979d0/
cc	data eps/1.0d-8/
cc	data eps/1.0d-0/
c	eps=xx
ccc	pip2=pi*0.5d0
c	pip4=pi*0.25d0
ccc	pip2p=pip2+eps
c	pip4p=pip4+eps
ccc	pip2=pi*1.5d0
c	pip4=pip4+pi
ccc	pip2m=pip2+eps
c	pip4m=pip4+eps
c	if(cdabs(e).lt.1.d-20) e=(1.0d-20,0.0d0)
c	x=dreal(e)
c	y=dimag(e)
c	fi=datan2(y,x)
c	if(y.lt.0.d0) fi=twopi+fi
cc
cc	fi is the arg(e) changing continuously from 0 to 2*pi
cc
c	psi=0.5d0*fi
cc
cc	psi is arg(k) changing continuously from 0 to pi
cc
ccc	if(psi.gt.pip2p.and.psi.lt.pip2m) psi=psi-pi
c	if(psi.gt.pip4p.and.psi.lt.pip4m) psi=psi-pi
c        if(psi.lt.psimin)psimin=psi
cc
cc	here psi is changing from 0 to pi/2+eps then jumps to
cc	-pi/2+eps and increases continuously till zero
cc	eps=1.d0
cc
c	r=dsqrt(cdabs(e))
c	x=dcos(psi)*r
c	y=dsin(psi)*r
c	cwn=dcmplx(x,y)
c	return
c	end	
cc040503        
C============SUBROUTINE EQMESH===================
        subroutine eqmesh(y,b,ypmin,ypmax,sol)
        implicit double complex (a-h,o-q,s-z)
        implicit double precision (r)
	parameter(maxs=250,nseg=20,nmax=2000,neqd=nmax)
	dimension y(0:*),b(1),sol(0:*)
	common /abscis/ rxeq(0:nmax),rx(0:neqd),nbo,numeq,numneq
        call splinc(rx(nbo),y(0),numneq-nbo,ypmin,ypmax,b)
	do 1 int=nbo+1,numeq
        call spintc(rx(nbo),y(0),b,numneq-nbo,rxeq(int),sol(int))
1	continue	
	return
	end
C============SUBROUTINE SCATWF===================
        subroutine scatwf(nbo,ck,match,sy,syp,dy1,dyp1,dy2,dyp2)
        implicit double complex(a-h,o-q,s-z)
        implicit double precision(r)        
        parameter(mni=3000,mnia=300)
        dimension cfc(0:3),cgc(0:3),cgcp(0:3),cfcp(0:3),sig(0:3)
	dimension sy(0:*),syp(0:*),dy1(0:*),dyp1(0:*),dy2(0:*),
     *  dyp2(0:*)
	dimension syk(mni),sykp(mni),x1(2)
        common/poti1/v1s,v10,v11,v12
        common/poti2/pas(mni),v(mni),v1(mni),v2(mni),intf
        common/genwork/par(3),rl
	common/radial/ renfac,retak,lorb,jorb
	rpi=4.d0*atan(1.d0)
	ci=(0.d0,1.d0)
c	write(*,*) 'nbo(soli1):',nbo
c        ener=ck**2/renfac
c	write(*,*) 'ck,e=',ck,ener
	rh=dreal(par(1))/nbo
       call soli1(rh,rl,sy(1),syp(1),syk,sykp,ck,v1s,
     * v10,v11,v12,nbo)    
        isens=1
c  integrate outwards and store in sy
	do 1 int=1,match
	int1=int+nbo
	sy(int1)=sy(int1-1)
	syp(int1)=syp(int1-1)
        call slix(isens,int,int,ck,pas,v,v1,v2,sy(int1),syp(int1),
     *  syk(int1),sykp(int1))
c	write(*,55) int1,sy(int1),syp(int1)
55	format(2x,'int=',i4,2x,4d12.4)
1	continue
	crho=ck*par(2)
	eta=retak/ck
	mode1=1
	kfn=0
	zlp=rl
	nlp=1
	ifail=1
	iprint=1
c
c	calculating coulomb functions at complex arguments
c
	call wclbes(crho,eta,zlp,nlp,cfc,cgc,cfcp,cgcp,sig,
     *  kfn,mode1,ifail,iprint)
c  O=G+i*F and O'
        out=cgc(0)+ci*cfc(0)
        outp=cgcp(0)+ci*cfcp(0)

c dF/dr & dO/dr
	cfcp(0)=ck*cfcp(0)
        outp=outp*ck

	intf1=nbo+intf
c integrate inward F(r) and store in dy1
	dy1(intf1)=cfc(0)
	dyp1(intf1)=cfcp(0)
	dyk=(0.d0,0.d0)
	dykp=(0.d0,0.d0)
	isens=-1
	do 2 int=intf,match+1,-1
	intf1=nbo+int
	dy1(intf1-1)=dy1(intf1)
	dyp1(intf1-1)=dyp1(intf1)
        if(intf.ge.match+1)call slix(isens,
     *int,int,ck,pas,v,v1,v2,dy1(intf1-1),dyp1(intf1-1),dyk,dykp)
2	continue
c dG/dr
c  integrate inward G(r) and store in dy2
	cgcp(0)=ck*cgcp(0)
	intf1=nbo+intf
	dy2(intf1)=cgc(0)
	dyp2(intf1)=cgcp(0)
	dyk=(0.d0,0.d0)
	dykp=(0.d0,0.d0)
	do 3 int=intf,match+1,-1
	intf1=nbo+int
	dy2(intf1-1)=dy2(intf1)
	dyp2(intf1-1)=dyp2(intf1)
        if(intf.ge.match+1)call slix(isens,
     *int,int,ck,pas,v,v1,v2,dy2(intf1-1),dyp2(intf1-1),dyk,dykp)
3	continue
	intf1=intf1-1

	x1(1)=(sy(intf1)*dyp2(intf1)/dy2(intf1)-syp(intf1))/
     +        (dy1(intf1)*dyp2(intf1)/dy2(intf1)-dyp1(intf1))
	x1(2)=(sy(intf1)*dyp1(intf1)/dy1(intf1)-syp(intf1))/
     +        (dy2(intf1)*dyp1(intf1)/dy1(intf1)-dyp2(intf1))
	a_u=0.5d0*(x1(2)+ci*x1(1))
	b_u=0.5d0*(x1(2)-ci*x1(1))
c below is the S-matrix of H_0
	smatrix=-b_u/a_u
	
c below is the corresponding phase shift of H_0
cc      after normaliaztion, u=sqrt(2/(pi S) ) i/2 (I-S O)
cc      u~a(I+b/a O)=x1 f + x2 g
	delta=(0.d0,-0.5d0)*log(smatrix)
	om=(0.d0,-1.d0)*delta
	tmp1=exp(om)
	ecm=ck*ck/renfac
	xn=sqrt(ck/(rpi*ecm))*(0.d0,0.5d0)
c	coef=xn/x1(1)*tmp1
	coef=xn/a_u*tmp1
c        write(8,*)' coef=',coef
c writes the S-matrix and delta for H_0
c	write(8,51) smatrix,delta
c	write(*,51) smatrix,delta

51	format('smatrix=',2d15.7,1x,'phaseshift=',2d15.7)
c	write(8,*) 'abs(smatrix) =',abs(smatrix)
c	pause
	nbmat=nbo+match
	nbintf=nbo+intf
	do 4 int=nbmat,nbintf
c	sy(int)=x1(1)*dy2(int)+x1(2)*dy1(int)
c	syp(int)=x1(1)*dyp2(int)+x1(2)*dyp1(int)
	sy(int)=x1(1)*dy1(int)+x1(2)*dy2(int)
	syp(int)=x1(1)*dyp1(int)+x1(2)*dyp2(int)
4	continue
	sy(0)=(0.d0,0.d0)
	do 5 int=1,nbintf
	sy(int)=sy(int)*coef
	syp(int)=syp(int)*coef
5	continue
	return
	end
cc =========== subroutine potgen =====================
	subroutine potgen(lbas,itz)
	parameter (mshp=652,nene=100)
	parameter (nstmax=mshp,mkval=100)
	implicit double precision(a-h,o-z)
	double complex epsi,depth,energy,ene(nene),g(nene),resp(nene),
     1 resid(nene),res,s,fk2,fk,yxs,yxn
	double complex pot,apot,spot,aspot
	dimension nzp(2)
	common/sabsa/rtol,rmax,rstep,nstep,ll,j2,nzp,nkval
c?
	common/potp/diff,bigr,depth,energy,c1,zet,
     1 bigrc,vso,etak,nsurf,iz,l,jj
c?
	common/potent/ vc(mshp),wc(mshp),vs(mshp),ws(mshp),
     1	vcb(mshp),vcnt(mshp)
	double complex potval(nstmax)
	common/radial/ renfac,retak,lorb,jorb
	double complex par
	common/genwork/par(3),rl
c	double complex yspn,bspn,cspn,dspn       
c	double complex ysps,bsps,csps,dsps
	double complex ypmin,ypmax
	double complex yspn,bspn
	double complex ysps,bsps
	double complex ysp,yex,dify
c	common/potspn/yspn(mshp),bspn(mshp),cspn(mshp),dspn(mshp),
c     *  nspn
	common/potspn/yspn(mshp),bspn(mshp),nspn
c	common/potsps/ysps(mshp),bsps(mshp),csps(mshp),dsps(mshp),
c     *  nsps
	common/potsps/ysps(mshp),bsps(mshp),nsps
	common /spcalc/ ispn,isps
c040415
c	common/c3pot/npar(2),at,zt,v_0(2),r_0,a,v_so(2)
	common/exc3pot/at,zt,v_0(2),r_0,a,v_so(2)
cc040415
        common/form/v0,rnuc,anuc
	common/form1/v00,vsoo,rnucc,anucc
	common /rads/ rxx(mshp),rhx(mshp)
	common/rolc1/ rolc1
	npind=(-1)**lbas
	if(abs(npind).ne.1) then
	 write(*,*)'Something is wrong with the parity. STOP' 
	 write(2,*)'Something is wrong with the parity. STOP' 
	 stop
	endif
	if(npind.eq.1) then
	 v0=v_0(1) 
	 vso=v_so(1)
	else
	 v0=v_0(2) 
	 vso=v_so(2)
	endif
c	v0=v_0
c	vso=v_so
	v00=v_0(1)
	vsoo=v_so(1)
	xpow=at**(1.d0/3.d0)
	rnuc=r_0*xpow
	rnucc=rnuc
	anuc=a
	anucc=anuc
	fact=1.d0
	if(itz.eq.0)then
	 ap=1.0082d0
	 iz=1
	 it0=(at-2.d0*zt)/2.d0
	 fact=-dsqrt(0.5d0*dfloat(it0))
	else
	 iz=nzp(itz)
	 if(iz.eq.1) then
	  ap=1.0089d0
	 else
	  ap=1.0075d0
	 endif
	endif
c	bigrc=1.19d0*xpow
	bigrc=rnuc
c	write(*,*) 'bigrc=',bigrc
	zet=zt
	if(iz.eq.1)zet=0.d0
	red2m=at*ap/(at+ap)
c	write(*,*)'at,ap=',at,ap
	c1=0.04783258d0*red2m*fact
c        c1=0.0482d0*0.8
	rolc1=c1
	renfac=c1
c	write(8,*) 'c1=',renfac
	if(itz.ne.0)then
	 if(iz.eq.1)then
c	  write(8,*)' this is a calculation for neutrons'
c	  write(*,*)' this is a calculation for neutrons'
	 else if(iz.eq.2)then
c	  write(8,*)' this is a calculation for protons'
c	  write(*,*)' this is a calculation for protons'
	 else
	  write(*,*)' nzp should be 1 or 2'
	  stop
	 endif
	endif
	etak=0.71993d0*zet*c1
c	write(8,*)'etak=',etak
	retak=etak
c	write(*,*)' retak=',retak
c?
	call pnpot(652,nstep,rstep,rmax,at)
c?
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	nsps=nstep+1
	nspn=nstep+1
	isps=0
	ispn=0
	do i=1,nsps
	 rxx(i)=0.05d0+i*rstep
c	 rhx(i)=rstep
	 ysps(i)=spot(dcmplx(rxx(i)))
	 yspn(i)=pot(dcmplx(rxx(i)))
	enddo
	ypmin=(ysps(2)-ysps(1))/rstep
	ypmax=(ysps(nsps)-ysps(nsps-1))/rstep
c	call spln(rxx,ysps,ypmin,ypmax,nsps,bsps,csps,dsps,rhx)
	call splinc(rxx,ysps,nsps,ypmin,ypmax,bsps)
	ypmin=(yspn(2)-yspn(1))/rstep
	ypmax=(yspn(nspn)-yspn(nspn-1))/rstep
c	call spln(rxx,yspn,ypmin,ypmax,nspn,bspn,cspn,dspn,rhx)
	call splinc(rxx,yspn,nspn,ypmin,ypmax,bspn)
	do i=1,nstep
	 rxt=i*rstep
	 isps=1
	 ysp=spot(dcmplx(rxt))
	 isps=0
	 yex=spot(dcmplx(rxt))
	 dify=yex-ysp
c	 write(*,123) rxt,yex,ysp,dify
123	 format('x=',f6.3,6d12.4)
	 ispn=1
	 ysp=pot(dcmplx(rxt))
	 ispn=0
	 yex=pot(dcmplx(rxt))
	 dify=yex-ysp
c	 write(*,124) yex,ysp,dify
124	 format(8x,6d12.4)
	enddo
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
102		format(1x,' pole at omega=',2f15.9,' residuum=',2e12.4)
100		format(1x,' l,2j=',2i5,' nzp=',i3,' h=',f10.4,
     1 ' rmax=',f10.4,' nstep=',i6)
101		format(' e-start=',6f10.5)
1653	format('  v_0,r_0,a,v_so',f6.2,2f5.2,f6.2)
	return
	end
c	call pnpot(652,nstep,rstep,rmax,at)
cc =========== subroutine pnpot ====================
	subroutine pnpot(n,m,h,rmx,at)
	parameter (mshp=652)
	implicit double precision(a-h,o-z)
	double complex energy,depth
	dimension vcb(n),vcnt(n)
	dimension vterm(mshp)
	common/form/v0,rnuc,anuc
	common/potp/diff,bigr,depth,energy,c1,zet,
     1 bigrc,vso,etak,nsurf,iz,l,jj
c***********************************************************************
c
c      subroutine to calculate state independent local potential terms
c
c
c      input parameters:
c
c      rc      coulomb radius
c      etak      coulomb parameter*wave number
c      h      step length in r
c      c1      2mu/hbar**2
c      m      number of points in r
c
c
c      results:
c
c***********************************************************************
      r=0.0
      nmin=1
      c=1.43986d0*zet
      bigr=rnuc
      diff=anuc
      v=v0
      depth=v
      nsurf=bigr/h+10
      return
      end
c**************************************************************************
C============FUNCTION SPOT=====================
	double complex function spot(x)
c  calculation of the singular potential which behaves as 1/r
	implicit double precision (a-h,o-z)
        parameter (mshp=652)
	double complex energy,bigk,depth,vvv,x
c	double complex ysps,bsps,csps,dsps
	double complex ysps,bsps
c	common/potsps/ysps(mshp),bsps(mshp),csps(mshp),dsps(mshp),
c     *  nsps
	common/potsps/ysps(mshp),bsps(mshp),nsps
	common /spcalc/ ispn,isps
	common/radial/ renfac,retak,lorb,jorb
	common /rads/ rxx(mshp),rhx(mshp)
	common/potp/diff,bigr,depth,energy,c1,zet,
     1  bigrc,vso,etak,nsurf,iz,ll,jj
        jj=jorb
        ll=lorb
c        write(*,*)'renfac,retak=',renfac,retak
c        write(*,*)'jj,ll=',jj,ll
	if(isps.eq.0) then
c
c   non-tabulated forms
c 
c  spin-orbit term
        if(jj.eq.0)then
        als=0.d0
        vls=0.d0
        else
	t=dreal(x)
	e1=exp((bigr-t)/diff)
	als=0.25d0*(jj*(jj+2)/4.d0-ll*(ll+1)-0.75d0)
	vls=als*4.d0/diff*e1/(1+e1)**2*vso
        endif
	vcou=0.d0
c  point Coulomb term outside
        if(iz.gt.1)then
	cost=1.43986d0*zet
	if(t.gt.bigrc)vcou=cost
        endif

	spot=c1*(vcou-vls)
	else
c  tabulated and splined term
	rxt=dreal(x)
        call spintc(rxx,ysps,bsps,nsps,rxt,spot)
	endif
	return
	end
***************************************************************************
C============FUNCTION POT=====================
	double complex function pot(x)
c   non-singular potential terms
	implicit double precision(a-h,o-z)
        parameter (mshp=652)
	double complex energy,bigk,depth,vvv,x
	common /potp/diff,bigr,depth,energy,c1,zet,
     1 bigrc,vso,etak,nsurf,iz,ll,jj
	double complex yspn,bspn
	common/potspn/yspn(mshp),bspn(mshp),nspn
	common /spcalc/ ispn,isps
	common /rads/ rxx(mshp),rhx(mshp)

	if(ispn.eq.0) then
c  non tabulated forms
c  W-S term
	t=dreal(x)
        if(abs(diff).lt.1.e-2)then
        vvv=depth
        if(t.gt.bigr)vvv=0.d0
        else
	e1=exp((bigr-t)/diff)
	vvv=depth*e1/(1.d0+e1)
        endif
c  Coulomb term
	vcou=0.d0
        if(iz.gt.1)then
	cost=1.43986d0*zet
	if(t.le.bigrc)vcou=cost/2.d0*(3.d0/bigrc-t*t/bigrc**3)
        endif

	pot=c1*(-vvv+vcou)
	else
c  tabulated and splined terms
	rxt=dreal(x)
        call spintc(rxx,yspn,bspn,nspn,rxt,pot)
	endif
	return
	end
***************************************************************************
C============FUNCTION ASPOT====================
	double complex function aspot(x)
c singular potential in the asymptotic region
c  Coulomb term
	implicit double precision (a-h,o-z)
        parameter (mshp=652)
	double complex energy,bigk,depth,vvv,x
	common/potp/diff,bigr,depth,energy,c1,zet,
     1  bigrc,vso,etak,nsurf,iz,ll,jj
	t=dreal(x)
	vcou=0.d0
        if(iz.gt.1)then
	cost=1.43986d0*zet
	if(t.gt.bigrc)vcou=cost
        endif
	aspot=c1*vcou
	return
	end
***************************************************************************
C============FUNCTION APOT=====================
	double complex function apot(x)
c  non-singular term in the asymptotic region: none
	implicit double precision(a-h,o-z)
	apot=(0.d0,0.d0)
	return
	end
*****************************************************************************
C============SUBROUTINE SLI======================
        subroutine sli(rtol,pas,v,v1,v2,par,rlorb,spot,pot,intf)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       this sbr. computes the partition and the associated vectors
c    (v,v1,v2) to be further used by sbr. slix which solves the
c    initial value problem for the radial schroedinger equation on i2
c    and on some subinterval of i3.
c       the equation to be solved is 
c                        y" = (v(x)- k**2)*y
c       where the complex function v(x) has the form
c                
c       v(x) = l*(l +1)/x**2 +spot(x)/x +pot(x)
c
c
c       input data(to be declared in the calling program):
c          rtol   accuracy requested in the results; rtol is a real double
c                 precision number
c          par    double complex vector with two components. par(1) and par(2)
c                 are the ends of the x domain of the system. the domain is
c                 supposed to be the segment in the complex plane which joints
c                 par(1) and par(2).
c          rlorb  double precision number for the orbital momentum l.
c
c       externals furnished by the user:
c          spot   function of the form
c                        double complex function spot(x)
c                 which furnishes spot(x), for any given double complex x.
c          pot    function of the form
c                               double complex function pot(x)
c                 which allows calculating  pot(x) for any given x
c
c       input/utput parameter:
c          intf   in input this furnishes the maximum number of intervals
c                 assigned to the partition. in output it furnishes the
c                 total number of intervals of the resultant partition 
c                 consistent with the required tolerance rtol. if
c                 it is found that this number exceeds the input value
c                 the program displayes the message
c                       "insufficient number of steps."
c          pas    double complex vector; pas(int) is the stepsize of the
c                 int-th interval of the partition.
c          v      double complex vector. value v(int), represents the gross
c                 part (conveniently normalized) of the potential on the 
c                 interval int.   
c          v1     double complex vector. value v1(int) represents the 
c                 (conveniently normalized) linear correction in the 
c                 potential on the interval int.
c          v2     the same as v1 for the quadratic correction.
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        implicit double complex (a-h,o-q,s-z)
        implicit double precision (r)
c        dimension par(1),pas(1),v(1),v1(1),v2(1)
        dimension par(2),pas(1),v(1),v1(1),v2(1)
	external spot,pot
        rlotol=rtol**(7.d0/6.d0)
        rl=(1.d-4)**(7.d0/6.d0)
        fac=(1.d2*rl/rlotol)**(1.d0/7.d0)
        nm=fac
        rnm=nm
        rld=abs(par(2)-par(1))
        versor=(par(2)-par(1))/rld
        int=0
        x1=par(1)
        pas(1)=versor
        ie=0
10      int=int+1
        icont=0
12      continue
	call sprel(x1,pas(int),v(int),v1(int),v2(int),spot,pot,rlorb)
        icont=icont+1
55      format(2x,'int=',i3,' icont=',i3,' x1=',2(f7.4),' pas=',
     *2(f7.4))    
        if(ie.eq.1) go to 15
        rv1=abs(v1(int)*v1(int))
        rhi=abs(pas(int))
        rv1=rv1/(rhi*rhi*rhi*rhi*rhi*rhi)
        rvbar=1.d-1
        rw=rv1*rvbar
        rdh=.1d0*rhi
        call hsize(rw,rdh,rh,rl)
        rdist=abs(x1+versor*rh-par(1))
        if(rdist.gt.rld-0.01*rh) ie=1
        if(int.eq.intf) ie=1
        if(int.eq.intf) write(*,*)' insufficient number of steps' 
        if(ie.eq.1) pas(int)=par(2)-x1
        if(ie.eq.1) go to 12
	if(icont.ge.3) go to 16
        if(abs(rh-rhi)/rhi.lt..11d0) go to 16
        pas(int)=versor*rh
        go to 12
16      x1=x1+pas(int)
        pas(int+1)=pas(int)
        go to 10
15      intf1=int
        ia=0
        if(abs(pas(intf1)).lt.1.d-12)ia=1
        if(ia.eq.1)intf1=intf1-1
        if(ia.eq.1)pas(intf1)=pas(intf1)+pas(intf1+1)
        intf2=intf1*nm
        if(intf2.gt.intf)write(*,*)'insufficient number of steps'
        if(intf2.gt.intf)stop
        intf=intf1
        do 1 int=1,intf
        inm=(int-1)*nm+1
        fac=pas(int)/rnm
        do 2 i=inm,inm+nm-1
        v(i)=fac
2       continue
1       continue
        intf=intf*nm
        do 3 int=1,intf
        pas(int)=v(int)
3       continue
        x1=par(1)
        do 20 k=1,intf
        call sprel(x1,pas(k),v(k),v1(k),v2(k),spot,pot,rlorb)
        x2=x1+pas(k)
        x1=x2
20      continue
        return
        end
**********************************************************************
C============SUBROUTINE HSIZE======================
        subroutine hsize(rw,rdh,rh,rlotol)
        implicit double precision (r)
        i=1
1       i=i+1
        if(i.ge.22) return
        m=0
        rh=rdh*(i-1)
        rh2=rh*rh
        rh6=rh2*rh2*rh2
        rw8=rw*rh6*rh2
        if(rw8/126.d1.le.rlotol) m=m+1
        if(rw8*rh/1134.d1.le.rlotol) m=m+1
        if(rw8/(9.d1*rh).le.rlotol) m=m+1
        if(m.eq.3) go to 1
        return
        end
********************************************************************
C============SUBROUTINE SPREL========================
          subroutine sprel(x1,h,v,v1,v2,spot,pot,rlorb)
          implicit double complex (a-h,o-q,s-z) 
          implicit double precision (r)
          dimension ppot(6),pspot(6)
	  external spot,pot
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c  this sbr. computes  v,v1,v2  on a single interval
c  whose ends are x1 and x1+h.          
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        rl=rlorb*(rlorb+1)
        call ssixval(pot,x1,h,ppot)
        call ssixval(spot,x1,h,pspot)
        do 3 m=1,6
        x=x1+.2d0*h*(m-1)
        ppot(m)=ppot(m)+pspot(m)/x+rl/(x*x)
    3 continue
        call sfitpot(ppot,h,v,v1,v2)
          return
          end
**************************************************************  
C============SUBROUTINE SSIXVAL=====================
        subroutine ssixval(func,x,h,pf)
        implicit double precision (r)
        implicit double complex (a-h,o-q,s-z)
        dimension pf(6)
        external func
        do 1 k=1,6
        pf(k)=func(x+.2d0*h*(k-1))
1       continue
        return
        end
***********************************************************************
C============SUBROUTINE SFITPOT===================
        subroutine sfitpot(v,h,v0,v1,v2)
        implicit double complex (a-h,o-q,s-z)
        dimension v(6)
        hh=h*h
      v0=hh*(19.d0*(v(1)+v(6))+75.d0*(v(2)+v(5))+5.d1*(v(3)+v(4)))/
     1288.d0
      v1=hh*(111.d0*(v(6)-v(1))+425.d0*(v(5)-v(2))
     1-15.d1*(v(4)-v(3)))/1344.d0
      v2=hh*25.d0*(5.d0*(v(1)+v(6))+6.d0*(v(2)+v(5))
     1-11.d0*(v(3)+v(4)))/1008.d0
      return
        end        
******************************************************************
C============SUBROUTINE SPI1======================
        subroutine spi1(x1,vc,v0,v1,v2)
        implicit double complex (a-h,o-q,s-z)
        implicit double precision (r)
        dimension cpfs(6),cpf(6)
        external spot,pot
        data czero/(0.d0,0.d0)/
        call ssixval(spot,czero,x1,cpfs)
        call sfitpot(cpfs,x1,vs0,vs1,vs2)
        call ssixval(pot,czero,x1,cpf)
        call sfitpot(cpf,x1,vn0,vn1,vn2)
        call copot(vs0,vs1,vs2,vn0,vn1,vn2,x1,vts,vt0,vt1,vt2)
        vc=vts/x1
        v0=vt0/(x1*x1)
        v1=vt1/(x1*x1*x1)
        v2=vt2/(x1*x1*x1*x1)
        return
        end
*********************************************************************
C============SUBROUTINE COPOT====================
        subroutine copot(vs0,vs1,vs2,v0,v1,v2,h,vts,vt0,vt1,vt2)
        implicit double complex (a-h,o-q,s-z)
        data two/(2.d0,0.d0)/,thr/(3.d0,0.d0)/,fou/(4.d0,0.d0)/,
     *  twe/(12.d0,0.d0)/
        vts=(vs0-two*(vs1-vs2))/h
        vt0=v0-two*(v1-v2)+fou*(vs1-thr*vs2)/h
        vt1=fou*(v1-thr*v2)+twe*vs2/h
        vt2=twe*v2
        return
        end

c *******************************************************************
c   subroutine wepc
c   purpose:
c     compute the limit of a sequence of complex numbers 
c     by the epsilon algorithm of p.wynn
c   method:
c     the epsilon algorithm of p.wynn
c   parameters:
c     input:
c     - s =    successive terms of the sequence to be accelerated
c     - inp =  number of calls made to the routine.
c              must be set to zero before first entry
c              and unchanged between calls.
c     - rtol = relative tolerance
c     - v    = working vector of dimension maxt
c     - maxt = the maximum number of terms of the sequence
c     output:
c     - cr=    the result (the rightmost term of the epsilon array)
c     - wr=    the precedent estimation of the limit (changed after
c              each successful calculation of the rightmost element
c              from the epsilon array belonging to an odd column) 
c     - rerr = the estimated relative error (a zero value means
c              that this estimation cannot be obtained)
c     - ier=   error flag
c       ier=0 - the tolerance was reached 
c       ier=1 - the tolerance was not reached in the frame of supplyed
c               terms
c       ier=2 - two successive terms are too close each other. the
c               calculation cannot continue. the last estimation of
c               the limit is returned as result.
c *******************************************************************
C============SUBROUTINE WEPC====================
        subroutine wepc(s,inp,rtol,v,maxt,cr,wr,rerr,ier)
        implicit double complex(a-h,o-q,s-z)
        implicit double precision(r)
        dimension v(maxt)
        data oner/(1.d0,0.d0)/,zero/(0.d0,0.d0)/,rsmall/1.d-7/
*********************************************************************
* reps is approximatively the smallest positive number of the machine
        data reps/1.0d-60/
*********************************************************************
        if(inp.ne.0) go to 1
        cr=s
        v(1)=s
        rerr=0.d0
        ier=1
        go to 5
1       k=0
        a1=s
        dif=a1-v(1)
        if(abs(dif).lt.reps) go to 8
        a0=oner/dif
        a2=a1
        a1=a0
        k=k+1
3       if(k.eq.inp) go to 2
        dif=a1-v(k+1)
        if(abs(dif).ge.reps) go to 7
8       ier=2
        rerr=0.d0
        go to 5
7       a0=v(k)+oner/dif
        v(k)=a2
        a2=a1
        a1=a0
        k=k+1
        go to 3
2       v(k)=a2
        v(k+1)=a1
        cr=v(k+1)
        ier=1
        if((inp/2)*2.eq.inp) go to 4
        if(inp.eq.1) go to 5
        rrvk=dreal(v(k))
        rivk=imag(v(k))
        rrvk2=dreal(v(k-2))
        rivk2=imag(v(k-2))
        rrtol=rtol
        if(abs(rrvk).gt.rsmall) rrtol=rtol*abs(rrvk)
        ritol=rtol
        if(abs(rivk).gt.rsmall) ritol=rtol*abs(rivk)
        if(abs(rrvk-rrvk2).le.rrtol.and.abs(rivk-rivk2).le.ritol) then
        ier=0
        cr=v(k)
        endif
        go to 5
4       rrcr=dreal(cr)
        ricr=imag(cr)
        rrvk1=dreal(v(k-1))
        rivk1=imag(v(k-1))
        rrtol=rtol
        if(abs(rrcr).gt.rsmall) rrtol=rtol*abs(rrcr)
        ritol=rtol
        if(abs(ricr).gt.rsmall) ritol=rtol*abs(ricr)
        if(abs(rrcr-rrvk1).le.rrtol.and.abs(ricr-rivk1).le.ritol) ier=0
5       inp=inp+1
        if(ier.eq.2) then
        cr=wr
        return
        endif
        if((inp/2)*2.eq.inp) return
        wr=cr
        return
        end

c   subroutine clevin(sofn,rofn,beta,n,arup,arlo,larray,estlim)
************************************************************************
*
*  description :
*  -------------
*
*  subroutine clevin computes levin's transformation.
*  (d. levin (j. comput. math., vol. b3, 1973, 371-388 ) using the
*  partial sums s(n) and terms a(n) of a given series as input data.
*
*  here, the definition of levin's sequence transformation according to
*  eq. (7.1-7) of e. j. weniger, comput. phys. rep., vol. 10 (1989),
*  189 - 371, is used.
*
*  in levin's sequence transformation it is assumed that the partial
*  sums s(n) of the series to be transformed can be written as
*
*  s(n) = s + r(n),
*
*  with s being the limit or antilimit of this series, and that the
*  exact remainder r(n) can be approximated by a remainder estimate
*  omega(n) in such a way that the ratio
*
*  [ s(n) - s ] / omega(n)
*
*  can be expressed as a by a poincare-type asymptotic expansion in
*  inverse powers of (n + beta).
*
*  levin suggested some simple remainder estimates omega(n) which can
*  be computed be determined from at most two terms a(n) of the series
*  to be transformed. with the help of these remainder estimates the
*  following variants of levin's sequence transformations result
*  (compare section 7.3 of e. j. weniger, comput. phys. rep.,
*  vol. 10 (1989), 189):
*
*  (n+beta) * a(n)                     : u transformation, eq. (7.3-5)
*  a(n)                                : t transformation, eq. (7.3-7)
*  a(n+1)                              : d transformation, eq. (7.3-7)
*  a(n + 1) * a(n) / [a(n) - a(n + 1)] : v transformation, eq. (7.3-11)
*
*  it is always assumed that the index n of the terms and the partial
*  sums s(n) satisfies n >= 0.
*
c  both numerator and denominator of this transformation are calculated 
c  with the help of a simple extension of a 2-dimensional nonlinear
c  3-term recurrence formula which was originally derived by 
c  t.fessler, w.f.ford and w.f.smith 
c  ( acm trans.math.software ,vol. 9 ( 1983 ), pp.346-354 ).           *
c                                                                      *
c  in this subroutine a variant of the so-called moving lozenge tech-  *
c  nique ( p.wynn, r.f.t.i.-chiffres, vol.9, pp.327-362 ( 1966 )) is   *
c  used, i.e., only one counterdiagonal of the numerator and denomina- *
c  tor tables have to be stored at a time. consequently, only two      *
c  1-dimensional arrays of appropriate length are needed in this       *
c  subroutine.                                                         *
c  in each call this subroutine calculates new counterdiagonals of the *
c  numerator and denominator tables and overwrites the previous entries.
c  in addition, a new estimate of the limit is calculated which is the *
c  in addition, a new estimate of the limit is calculated which is the *
c  ratio of the array elements arup(0) and arlo(0). hence, for every   *
c  new sequence element s-of-n this subroutine has to be called again.
c                                                                      *
c  input parameters :                                                  *
c  ----------------                                                    *
c                                                                      *
c  sofn   : element s(n) of the sequence to be accelerated             *
c                                                                      *
c  rofn   : estimate of the leading term of the poincare-type asympto- *
c           tic expansion of the remainder r(n) in inverse powers of 
c           (n+beta), i.e. omega(n).    
c           it is tacitly assumed that rofn is different from zero.    *
c                                                                      *
c  beta   : parameter contained in the inverse powers                  *
c           ( n + beta ) ** ( -j )                                     *
c           which describe the higher order contributions in the       *
c           poincare-type expansion of the remainder r(n).             *
c                                                                      *
c  n      : number of the last sequence element (n has to satisfy the  *
c           inequality n .ge. 0 )                                      *
c                                                                      *
c  arup   : 1-dimensional array to store the actual counterdiagonal of *
c           the numerator table.                                       *
c                                                                      *
c  arlo   : 1-dimensional array to store the actual counterdiagonal of *
c           the denominator table.                                     *
c                                                                      *
c  larray : length of the 1-dimensional arrays arup and arlo.          *
c                                                                      *
c  output parameter :                                                  *
c  ----------------                                                    *
c                                                                      *
c  estlim : approximation to the limit of the sequence s(n) which is   *
c           to be accelerated.                                         *
c           we have here estlim = arup(0)/arlo(0).                     *
c                                                                      *
c  machine-dependent parameters:                                       *
c  ----------------------------                                        *
c                                                                      *
c  huge,tiny : huge should be set close to but not identical with the  *
c              largest floating point number representable on the com- *
c              puter, and tiny should be set close to but not identical*
c              with the smallest floating point number representable   *
c              on the computer.                                        *
c                                                                      *
c***********************************************************************
C============SUBROUTINE CLEVIN==================
      subroutine clevin(sofn,rofn,beta,n,arup,arlo,larray,estlim)
      implicit double complex(a-h,o-z)
      double precision huge, tiny, one
      dimension arup(0:larray),arlo(0:larray)
      parameter ( huge = 1.d+60 , tiny = 1.d-60 , one = 1.d0 )
      zone=dcmplx(one,0.d0)
c
c  check whether n is less than zero
c
      if (n.lt.0) then
        write(*,1000)
1000    format('*** error in subroutine clevin ***'/
     1         '=================================')
        write(6,1010) n
1010    format('illegal input parameter n which must be .ge. zero'/
     1         'here we have n = ',i4)
        stop
c
c  check whether the arrays arlo and arup are large enough to store
c  all the elements required
c
      else if (n.gt.larray) then
        write(*,1000)
        write(*,1020) n,larray
1020    format('illegal input parameters n and larray'/
     1  'we must have n .le. larray'/
     2  'here we have n = ',i4,' and larray = ',i4)
        return
      end if
c
c  this is reached if the general levin transformation has to be
c  calculated.
c
c  first it has to be checked whether n .ge. 1 holds, i.e., whether we
c  proceed with the calculation of the next counter diagonals of the
c  two arrays or whether we have to initialize a new calculation.
c  this is done if n = 0 holds.
c
      arup(n)   = sofn / rofn
      arlo(n)   = zone  / rofn
c
c  check whether n .gt. 0 holds. in that case the recursive computation
c  of the numerators and denominators of the general levin transformatio
c  is to be continued.
c
      if ( n .gt. 0 ) then
        arup(n-1) = arup(n) - arup(n-1)
        arlo(n-1) = arlo(n) - arlo(n-1)
c
c  check whether n .gt. 1 holds. in that case the actual recursion
c  of the numerators and denominators starts.
c
        if (n.gt.1) then
          bn1     = beta + float(n-1)
          bn2     = beta + float(n)
          coef    = bn1 / bn2
          do 10 j = 2,n
            fact  = (beta+float(n-j)) * coef**(j-2) / bn2
c
c  perform the calculation of the new elements of the two tables.
c
            arup(n-j)  = arup(n-j+1) - fact*arup(n-j)
            arlo(n-j)  = arlo(n-j+1) - fact*arlo(n-j)
10        continue
        end if
      end if
c
c  compute the estimate for the limit of the sequence.
c
        if (abs(arlo(0)).lt.tiny) then
          estlim = dcmplx(huge,huge)
        else
          estlim = arup(0)/arlo(0)
        end if
      return
      end

c******************************************************************
C============SUBROUTINE ASSOL====================
	subroutine assol(x,pe,l,coul,u,up,uk,ukp)
c*****************************************************************
c     this sbr computes the outgoing solution o_l = g_l + i f_l  *
c     of the coulomb equation using an asymptotic series.        *
c     the summation of the asymptotic series is done by two      *
c     summation methods (wynn and levin) and the most accurate   *
c     solution is retained.                                      *
c     called routines: clevin,wepc
c
c     input parameters:                                          *
c     x = radius                                                 *
c     pe = k-value                                               *
c     l = angular momentum                                       *
c     coul = 2*k*eta, where eta is the sommerfeld parameter.     *
c     output parameters:                                         *
c     u = solution                                               *
c     up = d(u)/dr                                               *
c     uk = d(u)/dk                                               *
c     ukp= d(d(u)/dk)/dr                                         *
c*****************************************************************
      implicit double precision(r)
      implicit double complex(a-h,o-q,s-z)
	parameter (max=30,max1=max-1)
        dimension p(4,2,max),s(4,2),q(4,2),sq(4,2),qq(4,2),sqq(4,2)
     *  ,qw(4,2),sqw(4,2)
	dimension vvv(max)
      	dimension arup(0:max),arlo(0:max)
	data rtol/1.d-15/,reps/1.d-12/,rsmall/1.d-7/
	alpha=0.5d0*coul
 	eta=alpha/pe
        or=dcmplx(1.d0,0.d0)
        oi=dcmplx(0.d0,1.d0)
	oz=dcmplx(0.d0,0.d0)
	nmax=max
	nm1=nmax-1
	beta=or
        z=x*pe
        f1=alpha/(pe*z)
        f2=or-f1
	f3=oi*alpha/(pe*pe)
	f4=f3*pe
        f5=-.5d0*alpha/(pe*pe)
	zlog2=log(2.d0*z)
	s(1,1)=exp(oi*(z-eta*zlog2))
	s(2,1)=oi*f2*s(1,1)
	s(3,1)=(oi*f1/z-f2*f2)*s(1,1)
	s(4,1)=(-f1*(2.d0*oi/z+3.d0*f2)/z-oi*f2*f2*f2)*s(1,1)
	fs=f3*zlog2
	s(1,2)=fs*s(1,1)
	s(2,2)=(f3/z+oi*f2*fs)*s(1,1)
	s(3,2)=(-f1*(oi/z+2.d0*f2)/pe+(oi*f1/z-f2*f2)*fs)*s(1,1)
	w1=-.5d0*oi*(f4+l+or)*(f4-l)
	w2=f5*(2.d0*f4+or)
	p(1,1,1)=w1
	p(2,1,1)=-p(1,1,1)
	p(3,1,1)=2.d0*p(1,1,1)
	p(4,1,1)=-6.d0*p(1,1,1)
	p(1,2,1)=w2
	p(2,2,1)=-p(1,2,1)
        p(3,2,1)=2.d0*p(1,2,1)
	do 1 n=2,nmax
	w1=-.5d0*oi*(f4+l+n)*(f4-l+n-or)/float(n)
	w2=f5*(2.d0*(f4+n)-or)/float(n)
	p(1,1,n)=w1*p(1,1,n-1)
	p(2,1,n)=-n*p(1,1,n)
	p(3,1,n)=-(n+or)*p(2,1,n)
	p(4,1,n)=-(n+2.d0*or)*p(3,1,n)
	p(1,2,n)=w1*p(1,2,n-1)+w2*p(1,1,n-1)
	p(2,2,n)=-n*p(1,2,n)
	p(3,2,n)=-(n+or)*p(2,2,n)
1	continue
	zi=or/z
	zc=zi
	do 11 n=1,nmax
	do 10 i=1,4
	do 10 j=1,2
	if(i+j.gt.5) go to 10
	p(i,j,n)=p(i,j,n)*zc
10	continue
	zc=zc*zi
11	continue
	do 20 i=1,4
	do 20 j=1,2
	if(i+j.gt.5) go to 20
	inp=0
	ss=oz
	do 31 n=1,nm1
	term=p(i,j,n)
	ss=ss+term
31	continue
	sum=oz
	do 30 n=1,nm1
	term=p(i,j,n)
	sum=sum+term
	sofn=sum
	zofn=p(i,j,n+1)
	if(abs(zofn).le.rtol) go to 21
c  convergence acceleration by levin	
        call clevin(sofn,zofn,beta,inp,arup,arlo,max,estlim)
	if(inp.eq.0) then
	estlio=estlim
	go to 23
	endif
	rres=dreal(estlim)
	ries=imag(estlim)
	rreso=dreal(estlio)
	rieso=imag(estlio)
	rreps=reps
	if(abs(rres).gt.rsmall) rreps=reps*abs(rres)
	rieps=reps
	if(abs(ries).gt.rsmall) rieps=reps*abs(ries)
        if(abs(rres-rreso).le.rreps.and.abs(ries-rieso).le.rieps) 
     *  go to 21     
	estlio=estlim
c  epsilon algorithm by wynn
23	call wepc(sum,inp,reps,vvv,max,cr,wr,rerr,ier)
	if(ier.ne.1) go to 21
30 	continue
c	write(*,*) 'tol is not achieved'
21	continue 
	q(i,j)=wr
	qq(i,j)=ss
	qw(i,j)=estlim
20      continue
        zc=zi
	do 40 i=2,4
	do 41 j=1,2
	if(i+j.gt.5) go to 41
	q(i,j)=q(i,j)*zc
	qq(i,j)=qq(i,j)*zc
	qw(i,j)=qw(i,j)*zc
41	continue
	zc=zc*zi
40	continue
	q(1,1)=q(1,1)+or
	qq(1,1)=qq(1,1)+or
	qw(1,1)=qw(1,1)+or
	do 4 i=1,4
	do 5 j=1,2
	if(i+j.gt.5) go to 5
	sq(i,j)=dcmplx(0.d0,0.d0)
	sqq(i,j)=dcmplx(0.d0,0.d0)
	sqw(i,j)=dcmplx(0.d0,0.d0)
	do 6 k=1,i
	do 6 m=1,j
	if(k+m.gt.5) go to 6
	do 7 l1=1,i
	do 7 n=1,j
	if(l1+n.gt.5) go to 7
	if(k+l1-i.ne.1) go to 7
	if(m+n-j.ne.1) go to 7
	factor=or
   	if(i-k.eq.1.or.i-l1.eq.1) factor=(i-1)*or
        sq(i,j)=sq(i,j)+factor*s(k,m)*q(l1,n)
        sqq(i,j)=sqq(i,j)+factor*s(k,m)*qq(l1,n)
        sqw(i,j)=sqw(i,j)+factor*s(k,m)*qw(l1,n)
7	continue
6	continue
5	continue
4	continue
	u=sq(1,1)
	up=pe*sq(2,1)
        upp=pe*pe*sq(3,1)
	v=sqq(1,1)
	vp=pe*sqq(2,1)
        vpp=pe*pe*sqq(3,1)
	w=sqw(1,1)
	wp=pe*sqw(2,1)
        wpp=pe*pe*sqw(3,1)
	uk=sq(1,2)+x*sq(2,1)
	ukp=sq(2,1)+pe*sq(2,2)+z*sq(3,1)
c	ukpp=pe*(2.d0*sq(3,1)+pe*sq(3,2)+z*sq(4,1))
	vk=sqq(1,2)+x*sqq(2,1)
	vkp=sqq(2,1)+pe*sqq(2,2)+z*sqq(3,1)
c	vkpp=pe*(2.d0*sqq(3,1)+pe*sqq(3,2)+z*sqq(4,1))
	wk=sqw(1,2)+x*sqw(2,1)
	wkp=sqw(2,1)+pe*sqw(2,2)+z*sqw(3,1)
c	wkpp=pe*(2.d0*sqw(3,1)+pe*sqw(3,2)+z*sqw(4,1))
	pot=(l*(l+1.d0)/x+2.d0*eta*pe)/x-pe*pe
	dpotk=-2.d0*pe
	du=upp-pot*u
	dv=vpp-pot*v
	dw=wpp-pot*w
	rdv=abs(dv)
	rdu=abs(du)
	rdw=abs(dw)
	rdmin=min(rdv,rdu,rdw)
	if(rdmin.eq.rdu) go to 25
	if(rdmin.eq.rdv) then
	u=v
	up=vp
	uk=vk
	ukp=vkp
	else
	u=w
	up=wp
	uk=wk
	ukp=wkp
	endif
25	quady2=.5d0*(ukp*u-up*uk)/pe
c
c	quad2 is the integral from infinity to x of u*u; u is the
c	outgoing coulomb function
c
	return
	end
********************************************************************
C============SUBROUTINE SOLI3=========================        
        subroutine soli3(ck,dy,dyp,dyk,dykp)
c      
c    this computes the outgoing solution (dy),its derivative with respect to
c    x (dy),the derivative of y against k (dyk) and the derivative of dy
c    against k (dykp) at point par(3), the l.e.end of interval i2.
c      
c       input parameter
c       ck   double complex value for momentum k
c
c       other data are communicated through 'common' blocks
c
        implicit double complex(a-h,o-q,s-z)
        implicit double precision(r)
        parameter(mni=3000,mnia=300)
        common/poti3/apas(mnia),av(mnia),av1(mnia),av2(mnia),intaf
        common/genwork/par(3),rl
        isens=-1
        coul=aspot(par(3))
        l=rl
        call assol(par(3),ck,l,coul,dy,dyp,dyk,dykp)
        call slix(isens,intaf,1,ck,apas,av,av1,av2,dy,dyp,dyk,dykp)
        return
        end
C============SUBROUTINE SFINDEV======================        
        subroutine sfindev(ck,match,swght,gap,dwght,rtol,ierr)
c
c     this computes the eigenvalue of the boundary value problem. newton
c     method is used in the iteration
c
c       input parameter
c       match   the interval at the r.h.end of which the solutions
c               should be matched
c
c       input/output parameter:
c       ck      in input it furnishes a guess value. in output, it contains
c               the resultant eigenvalue.
c
c       output parameters:
c       swght   double complex numbers which are further on used for the
c       gap     normalization of the eigenfunction
c       dgwht
c           
        implicit double complex(a-h,o-q,s-z)
        implicit double precision(r)        
        parameter(mni=3000,mnia=300)
	parameter(maxnewit=20)
        common/poti1/v1s,v10,v11,v12
        common/poti2/pas(mni),v(mni),v1(mni),v2(mni),intf
        common/poti3/apas(mnia),av(mnia),av1(mnia),av2(mnia),intaf
        common/genwork/par(3),rl
	dimension y(1),yp(1),yk(1),ykp(1)
	
        newit=0
1       call soli1(dreal(par(1)),rl,y,yp,yk,ykp,ck,v1s,v10,v11,
     *v12,1)    
	sy=y(1)
	syp=yp(1)
	syk=yk(1)
	sykp=ykp(1)
        isens=1
        call slix(isens,1,match,ck,pas,v,v1,v2,sy,syp,syk,sykp)
        call soli3(ck,dy,dyp,dyk,dykp)
        isens=-1
        if(intf.ge.match+1)call slix(isens,
     *intf,match+1,ck,pas,v,v1,v2,dy,dyp,dyk,dykp)
        phi=sy*dyp-dy*syp
        dphi=syk*dyp+sy*dykp-dyk*syp-dy*sykp
        dphi=(dphi*sy*dy-phi*(syk*dy+sy*dyk))/(sy*dy)
        dk=-phi/dphi
c        write(*,*)' increment dk is ',dk
        ck=ck+dk
c        write(*,*)' new k is ',ck
        newit=newit+1
	if(abs(dk/ck).le.rtol) then
	ierr=0
	go to 2
	else
        if(newit.le.maxnewit) go to 1
	ierr=1
	go to 2
	endif
2       swght=.5d0*(syp*syk-sykp*sy)/ck
        dwght=.5d0*(dyp*dyk-dykp*dy)/ck
        gap=sy/dy
        return
        end
C============SUBROUTINE SNORM=======================
        subroutine snorm(nbo,ck,gap,match,swght,dwght,y,yp)
c
c       this computes the normalized eigenfunction associated with the
c       eigenvalue ck
c
c       input parameters:
c       nbo    number of equidistant points on i1 at which the solution
c              is needed
c       ck     the eigenvalue, as resulted from sbr.sfindev
c       gap    three double complex numbers as resulted from sbr. sfindev
c       swght
c       dwght        
c
c 
c        output parameters
c        y      double complex vector
c        
        
        implicit double complex(a-h,o-q,s-z)
        implicit double precision(r)        
        parameter(mni=3000,mnia=300)
        common/poti1/v1s,v10,v11,v12
        common/poti2/pas(mni),v(mni),v1(mni),v2(mni),intf
        common/poti3/apas(mnia),av(mnia),av1(mnia),av2(mnia),intaf
        common/genwork/par(3),rl
        dimension y(0:*),yp(0:*)
        dimension yk(10),ykp(10)
        y(0)=0.d0
        qsd=swght-dwght*gap*gap
        sqsd=sqrt(qsd)
51      format(2x,'sw=',2d16.8,2x,'dw=',2d16.8)
        swght=(1.d0,0.d0)/sqsd  
        dwght=gap/sqsd
c        write(8,51) swght,dwght
	rh=dreal(par(1))/nbo
        call soli1(rh,rl,y(1),yp(1),yk,ykp,ck,
     *v1s,v10,v11,v12,nbo)
        do 1 n=1,nbo
        y(n)=y(n)*swght
        yp(n)=yp(n)*swght
        yk(n)=yk(n)*swght
        ykp(n)=ykp(n)*swght
1       continue
        isens=1
        do 25 int=1,match
        int1=int+nbo
        y(int1)=y(int1-1)
        yp(int1)=yp(int1-1)
        call slix(isens,int,int,ck,pas,v,v1,v2,y(int1),yp(int1),
     *yk(nbo),ykp(nbo))
25      continue
        call soli3(ck,dy,dyp,dyk,dykp)
        isens=-1
        int1=intf+nbo
        y(int1)=dy*dwght
        yp(int1)=dyp*dwght
        dyk=dyk*dwght
        dykp=dykp*dwght
        do 28 int=intf,match+1,isens
        int1=int+nbo
        y(int1-1)=y(int1)
        yp(int1-1)=yp(int1)
        call slix(isens,int,int,ck,pas,v,v1,v2,y(int1-1),
     *yp(int1-1),dyk,dykp)
28      continue
        return
        end
C============SUBROUTINE SOLI1=====================
      subroutine soli1(rh,rl,u,up,ue,uep,pe,vc,v0,v1,v2,nbo)
c
c       this computes the regular solution at a series of points on i1.
c
c       input parameters
c       rl          double precision number for the angular momentum l
c       vc,v0,v1,v2 double complex numbers as resulted from sbr.spi1    
c       pe          double complex number representing the momentum k
c       rh          double precision number for the r.h.end of i1      
c       nbo         number of equidistant points on i1 at which the solution
c                   is required.
c
c       output parameters
c       u           double complex vector. in u(i) the value of the solution
c                   at point x_i=i*rh is found. i ranges from 1 to nbo
c       up          the same for the first derivative of u
c       ue          the same for the first derivative of u against k
c       uep         the same for the first derivative of up against k
c
      implicit double precision(r)
      implicit double complex(a-h,o-q,s-z)
      dimension u(1),up(1),ue(1),uep(1)
      dimension f(4),ca(65),cap(65),a(65),ap(65),cfu(65),cfue(65)
      r=nbo*rh
      e=pe*pe
      l=rl
      rlmem=rl
      f(1)=r*vc
      f(2)=(v0-pe*pe)*r*r
      fp=-2.d0*pe*r*r
      f(3)=v1*r*r*r
      f(4)=v2*r*r*r*r
      cz=dcmplx(0.d0,0.d0)
      cu=dcmplx(1.d0,0.d0)
      iqd=0
      iqu=0
      iql=0
      do 20 i=1,65
      cfu(i)=cz
      cfue(i)=cz
      a(i)=cz
20    ap(i)=cz
      a(2)=cu
      cfu(2)=cu
c       write(8,1001) vc,pe
1001    format(2x,'vc=',2g14.6,2x,'k=',2g14.6)
      cpot=l*(l+1)/(r*r)+vc/r+(v0-pe*pe)+v1*r+v2*r*r
      do 10 kpert=1,15
      do 1 iq=iqd-1,iqu+2
      iqf=iq+3
      ca(iqf)=cz
      cap(iqf)=cz
      ipu=2
      ipl=iq-iqu
      if(ipl.le.-1) ipl=-1
      ipl=-1
      ipu=2
      if(ipu.ge.iq) ipu=iq
      do 2 ip=ipl,ipu
      ca(iqf)=ca(iqf)+f(ip+2)*a(iq-ip+2)
      cap(iqf)=cap(iqf)+f(ip+2)*ap(iq-ip+2)
2     continue
      cap(iqf)=cap(iqf)+fp*a(iq+2)
1     continue
      do 3 iq=1,iqu+4
      ap(iq)=cz
3     a(iq)=cz
      iql=iql+1
      iqu=iqu+4
      do 4 iq=iql,iqu
      den=(iq+l+1)*(iq+l)-l*(l+1)
      a(iq+2)=ca(iq+1)/den
      ap(iq+2)=cap(iq+1)/den
      cfu(iq+2)=cfu(iq+2)+a(iq+2)
      cfue(iq+2)=cfue(iq+2)+ap(iq+2)               
4     continue
      y=cu
      yp=(l+1)*cu
      yp=yp/r
      ye=cz
      yep=cz
      ypp=(l+1)*l*cu
      ypp=ypp/(r*r)
      yepp=cz
      do 11 iq=1,iqu
      iqc=iq+2
      su=(iq+l+1)*cfu(iqc)/r
      sue=(iq+l+1)*cfue(iqc)/r
      y=y+cfu(iqc)
      yp=yp+su
      ypp=ypp+(iq+l)*su/r
      ye=ye+cfue(iqc)
      yep=yep+sue
      yepp=yepp+(iq+l)*sue/r
11    continue
      rtol=1.d-15
50    format(2x,i4,6(d10.3,1x))
      crty=cpot*y
      crtye=-2.d0*pe*y+cpot*ye
      dy=ypp-crty
      dye=yepp-crtye
      if(abs(dy).le.rtol*abs(ypp).and.abs(dye).le.rtol*abs
     1(yepp)) go to 15
10    continue
c      write(*,*)' warning: the expected precision was not reached.'
c      pause
15    continue            
52    format(6(2x,d12.5))      
51    format(6x,12(d10.3,1x))
      do 16 n=1,nbo
        us=cz
        ups=cz
        upps=cz
        ues=cz
        ueps=cz
        uepps=cz
      rl=n*rh/r
      do 19 in=1,iqu
        iq=iqu-in+1
      iqc=iq+2
      su=(iq+l+1)*cfu(iqc)
      sue=(iq+l+1)*cfue(iqc)
      us=rl*(cfu(iqc)+us)
      ups=rl*(su+ups)
      upps=rl*((iq+l)*su+upps)
      ues=rl*(cfue(iqc)+ues)
      ueps=rl*(sue+ueps)
      uepps=rl*((iq+l)*sue+uepps)
19    continue
        u(n)=cu+us
        up(n)=(l+1)*cu+ups
        upp=l*(l+1)*cu+upps
        ue(n)=ues
        uep(n)=ueps
        uepp=uepps
      ril=1.d0
      do 21 i=1,l+1
      ril=ril*rl
21    continue
      rl=n*rh
      u(n)=u(n)*ril
      up(n)=up(n)*ril/rl
      upp=upp*ril/(rl*rl)
      ue(n)=ue(n)*ril
      uep(n)=uep(n)*ril/rl
      uepp=uepp*ril/(rl*rl)
      cpot=l*(l+1)/(rl*rl)+vc/rl+v0-e+v1*rl+v2*rl*rl
      crty=u(n)*cpot
      dy=upp-crty
      crtye=-2.d0*pe*u(n)+cpot*ue(n)
      dye=uepp-crtye
c      write(*,50)n,rl
c      write(*,51)u(n),ue(n)
c      write(*,51)up(n),uep(n)
c      write(*,51)upp,uepp
c      write(*,51)dy,dye
c      write(*,53)
53    format(/)      
16    continue                  
      rl=rlmem
      return
      end
C============SUBROUTINE SLIX=====================        
        subroutine slix(isens,inti,intf,ck,pas,v,v1,v2,y,yp,yk,ykp)
c
c       this sbr. propagates the solution along a prescribed interval
c
c       input parameters
c       isens     sense of the propagation; isens=1 for forward while
c                 isens=-1 for backward propagation.
c       inti      the label of the interval from which the propagation begins.
c       intf      the label of the last interval to be covered
c                 note: inti<intf requires isens=1,and so the solution is 
c                 propagated forwards from the l.h.end of inti up to the 
c                 r.h.end of intf;
c                 inti>intf requires isens=-1; the solution is then propagated
c                 backwards from the r.h.end of inti up to the l.h.end of intf.
c                 one can also consider the propagation on a single interval.
c                 then inti=intf and, if isens is settled 1, the propagation
c                 is forwards, from the l.h.end to the r.h.end of inti;
c                 if isens=-1, the propagation is backwards, from the r.h.end
c                 to the l.h.end of inti.
c       ck        the double complex value of the momentum k
c    pas,v,v1,v2  double precision vectors with the values resulted from
c                 sbr sli
c
c       input/output parameters
c   y,yp,yk,ykp   four double complex parameters. in input, they furnish the
c                 initial values for the solution (y), its first derivative
c                 against x (yp), its first derivative against k (yk)
c                 and the first derivative of yp against k (ykp).
c                 in output these parameters contain the value of the same
c                 quantities at the other end of the domain.
c
                  
        implicit double precision(r)
        implicit double complex(a-h,o-q,s-z)
        dimension pas(1),v(1),v1(1),v2(1)
        e=ck*ck
        pe=ck
        y0=y
        yp0=yp
        ye0=yk
        yep0=ykp
       do 2 n=inti,intf,isens
       h=pas(n)
       eh=e*h*h
      call tmcpm2(eh,v(n),v1(n),v2(n),csi,eta0,eta1,eta2,eta3,
     1eta4,eta5,eta6,t11,t12,t21,t22)
      if(isens.lt.0) go to 4
      y1=t11*y0+t12*yp0*h
      yp1=t21*y0/h+t22*yp0
      ye1=t11*ye0+t12*yep0*h-pe*h*h*(y0*eta0+yp0*h*eta1)
      yep1=t21*ye0/h+t22*yep0-pe*h*(y0*(csi+eta0)+yp0*h*eta0)
      go to 5
4     y1=t22*y0-t12*yp0*h
      yp1=-t21*y0/h+t11*yp0   
      ye1=t22*ye0-t12*yep0*h-pe*h*h*(y0*eta0-yp0*h*eta1)
      yep1=-t21*ye0/h+t11*yep0+pe*h*(y0*(csi+eta0)-yp0*h*eta0)
5     y0=y1
      yp0=yp1
      ye0=ye1
      yep0=yep1
2     continue
        y=y0
        yp=yp0
        yk=ye0
        ykp=yep0
        return
        end

C============SUBROUTINE TMCPM2================
      subroutine tmcpm2(eh,v,v1,v2,csi,eta0,eta1,eta2,eta3,
     1eta4,eta5,eta6,t11,t12,t21,t22)
c
c       calculation of the elements t11,..,t22 of the transfer
c         matrix t which is used to
c       advance the solution and its derivative from a point
c       r to r+h i.e. outward
c
c     (   y(r+h)  )      ( t11     t12*h  )   ( y(r) )
c     (           )   =  (                ) * (      )
c     (   y'(r+h) )      ( t21/h   t22    )   ( y'(r))
c
c       or from r to r-h i.e. inward
c
c     (   y(r-h)  )      ( t22     -t12*h  )   ( y(r) )
c     (           )   =  (                 ) * (      )
c     (   y'(r-h) )      ( t21/h   t11     )   ( y'(r))
c
      implicit double precision(r)
      implicit double complex(a-h,o-q,s-z)
        z=v-eh
        call cgebas(z,csi,eta0,eta1,eta2,eta3,eta4,eta5,eta6)
c
c       corrections of the 1-st order perturbation theory
c
        corr1=v1*eta1
        corr2=v2*eta2
        v11=v1*v1/6.d0
        v12=2.d0*v1*v2
        v22=.3d0*v2*v2
        auxv=v11+v22/3.d0
        auxc=auxv*eta2
c
c       corrections of the 2-nd order perturbation theory
c
        corr11=-auxc+(v12-v22)*eta3
        corr21=-auxv*eta1-(7.d0*v11+2.d0*v22)*eta2+19.d0*v22*eta3
        corr12=-auxv*eta3+3.d0*v22*eta4
        corr22=-auxc-(v12+v22)*eta3
c
c       calculation of the transfer matrix by adding corrections 
c       of the order 0,1 and 2 together
c
        t11=csi-corr1+corr11
        t21=z*(eta0+corr2)+corr21
        t12=eta0-corr2+corr12
        t22=csi+corr1+corr22
        return
        end
C============SUBROUTINE CGEBAS================= 
      subroutine cgebas(z,csi,eta0,eta1,eta2,eta3,eta4,eta5,eta6)
c
c      subroutine to calculate the quantities csi,..,eta6 defined in p.198
c      of ref.1.  by eq. (3.4.15)-(3.4.18b)
c
      implicit double precision(r)
      implicit double complex(a-h,o-q,s-z)
      rtz=1.d0
      ru=abs(z)
      if(ru.lt.rtz) go to 72
      squ=sqrt(z)
      squ=squ*dcmplx(0.d0,1.d0)
      eta0=sin(squ)/squ
      csi=cos(squ)
   71 eta1=(csi-eta0)/z
      eta2=(eta0-3.d0*eta1)/z
      eta3=(eta1-5.d0*eta2)/z
      eta4=(eta2-7.d0*eta3)/z
      eta5=(eta3-9.d0*eta4)/z
      eta6=(eta4-11.d0*eta5)/z
      return
   72 eta6=(269583552.d4+z*(89861184.d3+z*(1321488.d3+z*(11592.d3+
     1z*(69.d3+z*(3.d2+z))))))/3643017329952.d5
      eta5=(14018344704.d2+z*(539167104.d2+z*(89861184.d1+z*(
     1880992.d1+z*(5796.d1+z*(276.d0+z))))))/14572069319808.d3
      eta4=z*eta6+11.d0*eta5
      eta3=z*eta5+9.d0*eta4
      eta2=z*eta4+7.d0*eta3
      eta1=z*eta3+5.d0*eta2
      eta0=z*eta2+3.d0*eta1
      csi=z*eta1+eta0
      return
      end
c***************************************************************
C============SUBROUTINE SPLINC=================
      subroutine splinc(rx,y,n,yp1,ypn,y2)
c
c  calculation of the spline coefficients for knots points x(i),y(i),i=1,..,n
c  and store them in y2(1),..,y2(n)      
c
        parameter(mni=3000,mnia=300)
      implicit double complex (a-h,o-q,s-z)
      implicit double precision (r)
      dimension rx(n),y(n),y2(n),u(mni)
      if (dreal(yp1).gt..99e30) then
        y2(1)=0.d0
        u(1)=0.d0
      else
        y2(1)=-0.5d0
        u(1)=(3.d0/(rx(2)-rx(1)))*((y(2)-y(1))/(rx(2)-rx(1))-yp1)
      endif
      do 11 i=2,n-1
        rsig=(rx(i)-rx(i-1))/(rx(i+1)-rx(i-1))
        p=rsig*y2(i-1)+2.d0
        y2(i)=(rsig-1.d0)/p
        u(i)=(6.d0*((y(i+1)-y(i))/(rx(i+1)-rx(i))-(y(i)-y(i-1))
     *      /(rx(i)-rx(i-1)))/(rx(i+1)-rx(i-1))-rsig*u(i-1))/p
11    continue
      if (dreal(ypn).gt..99e30) then
        qn=0.d0
        un=0.d0
      else
        qn=0.5d0
        un=(3.d0/(rx(n)-rx(n-1)))*(ypn-(y(n)-y(n-1))/(rx(n)-rx(n-1)))
      endif
      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.d0)
      do 12 k=n-1,1,-1
        y2(k)=y2(k)*y2(k+1)+u(k)
12    continue
      return
      end
c******************************************************************	     
C============SUBROUTINE SPINTC===================
      subroutine spintc(rxa,ya,y2a,n,rx,y)
c
c  spline interpolation using the spline coefficients in y2a      
c
      implicit double complex (a-h,o-q,s-z)
      implicit double precision (r)
      dimension rxa(n),ya(n),y2a(n)
      klo=1
      khi=n
1     if (khi-klo.gt.1) then
        k=(khi+klo)/2
        if(rxa(k).gt.rx)then
          khi=k
        else
          klo=k
        endif
      goto 1
      endif
      rh=rxa(khi)-rxa(klo)
      if (rh.eq.0.d0) pause 'bad xa input.'
      ra=(rxa(khi)-rx)/rh
      rb=(rx-rxa(klo))/rh
      y=ra*ya(klo)+rb*ya(khi)+
     *((ra**3-ra)*y2a(klo)+(rb**3-rb)*y2a(khi))*(rh**2)/6.d0
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccc  FINALIZA MODULOS PARA CALCULAR LAS ccccccccccccccccccccccc
ccccccccccccccc   FUNCIONES DE SCATTERING Y GAMOW   ccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cc=======================================================================
cc======  FINISH THE PROGRAM WHICH SOLVE THE SCHRODINGER EQUATION =======
cc======                   FOR COMPLEX ENERGY                     =======
cc=======================================================================



c$$$	rphi=datan(1.d0)*2.d0-datan2(dimag(ck),dreal(ck))
c$$$	rk=abs(ck)
c$$$	rinf=25.d0/rk
c$$$	ci=cmplx(0.d0,1.d0)
c$$$	coul=2.d0*retak
c$$$	call gauleg(0,rinf,rgaussx,rgaussw,npoints)
c$$$	csum2=0.d0
c$$$	do i=1,npoints
c$$$	    cr=par(2)+rgaussx(i)*exp(ci*rphi)
c$$$	    call assol(cr,ck,lorb,coul,u,up,uk,ukp)
c$$$	    csum2=csum2+rgaussw(i)*up*up*dwght*dwght*exp(ci*rphi)
c$$$	end do
c$$$	write(*,*) csum1+csum2

	subroutine spline (x, y, b, c, d, n)
c!======================================================================
c!       Calculate the coefficients b(i), c(i), and d(i), i=1,2,...,n
c!       for cubic spline interpolation
c!       s(x) = y(i) + b(i)*(x-x(i)) + c(i)*(x-x(i))**2 + d(i)*(x-x(i))**3
c!       for  x(i) <= x <= x(i+1)
c!       Alex G: January 2010
c!----------------------------------------------------------------------
c!       input..
c!       x = the arrays of data abscissas (in strictly increasing order)
c!       y = the arrays of data ordinates
c!       n = size of the arrays xi() and yi() (n>=2)
c!       output..
c!       b, c, d  = arrays of spline coefficients
c!       comments ...
c!       spline.f90 program is based on fortran version of program spline.f
c!       the accompanying function fspline can be used for interpolation
c!======================================================================
	implicit none
	integer n
        double precision x(n)
	double complex y(n), b(n), c(n), d(n)
	integer i, j, gap
        double precision h

	gap = n-1
c       check input
	if ( n < 2 ) return
	if ( n < 3 ) then
	   b(1) = (y(2)-y(1))/(x(2)-x(1)) ! linear interpolation
	   c(1) = 0.
	   d(1) = 0.
	   b(2) = b(1)
	   c(2) = 0.
	   d(2) = 0.
	   return
	end if
c       
c       step 1: preparation
c       
        d(1) = x(2) - x(1)
        c(2) = (y(2) - y(1))/d(1)
        do i = 2, gap
	d(i) = x(i+1) - x(i)
	b(i) = 2.0*(d(i-1) + d(i))
	c(i+1) = (y(i+1) - y(i))/d(i)
	c(i) = c(i+1) - c(i)
	end do
c       
c       step 2: end conditions 
c       
	b(1) = -d(1)
	b(n) = -d(n-1)
        c(1) = 0.0
        c(n) = 0.0
	if(n /= 3) then
	   c(1) = c(3)/(x(4)-x(2)) - c(2)/(x(3)-x(1))
	   c(n) = c(n-1)/(x(n)-x(n-2)) - c(n-2)/(x(n-1)-x(n-3))
	   c(1) = c(1)*d(1)**2/(x(4)-x(1))
	   c(n) = -c(n)*d(n-1)**2/(x(n)-x(n-3))
	end if
c       
c       step 3: forward elimination 
c       
        do i = 2, n
	h = d(i-1)/b(i-1)
	b(i) = b(i) - h*d(i-1)
	c(i) = c(i) - h*c(i-1)
	end do
c       
c       step 4: back substitution
c       
        c(n) = c(n)/b(n)
        do j = 1, gap
	i = n-j
	c(i) = (c(i) - d(i)*c(i+1))/b(i)
	end do
c       
c       step 5: compute spline coefficients
c       
	b(n) = (y(n) - y(gap))/d(gap) + d(gap)*(c(gap) + 2.0*c(n))
        do i = 1, gap
	b(i) = (y(i+1) - y(i))/d(i) - d(i)*(c(i+1) + 2.0*c(i))
	d(i) = (c(i+1) - c(i))/d(i)
	c(i) = 3.*c(i)
	end do
        c(n) = 3.0*c(n)
        d(n) = d(n-1)
	end subroutine spline

	function cispline(u, x, y, b, c, d, n)
c!======================================================================
c!       function ispline evaluates the cubic spline interpolation at point z
c!       ispline = y(i)+b(i)*(u-x(i))+c(i)*(u-x(i))**2+d(i)*(u-x(i))**3
c!       where  x(i) <= u <= x(i+1)
c!----------------------------------------------------------------------
c!       input..
c!       u       = the abscissa at which the spline is to be evaluated
c!       x, y    = the arrays of given data points
c!       b, c, d = arrays of spline coefficients computed by spline
c!       n       = the number of data points
c!       output:
c!       ispline = interpolated value at point u
c!=======================================================================
	implicit none
        double complex cispline
	integer n
        double precision  u, x(n)
	double complex y(n), b(n), c(n), d(n)
	integer i, j, k
        double precision dx

c!       if u is ouside the x() interval take a boundary value (left or right)
	if(u <= x(1)) then
	   dx=u-x(1)
	   cispline = y(1)+ dx*(b(1) + dx*(c(1) + dx*d(1)))
	   return
	end if
	if(u >= x(n)) then
	   dx=u-x(n)
	   cispline = y(n)+ dx*(b(n) + dx*(c(n) + dx*d(n)))
	   return
	end if

c       *
c       binary search for for i, such that x(i) <= u <= x(i+1)
c       *
	i = 1
	j = n+1
        do while (j > i+1)
	k = (i+j)/2
	if(u < x(k)) then
	   j=k
	else
	   i=k
	end if
	end do
c       *
c       evaluate spline interpolation
c       *
        dx = u - x(i)
	cispline = y(i) + dx*(b(i) + dx*(c(i) + dx*d(i)))
	end function cispline
